Inelastic effects in molecular junctions in the Coulomb and Kondo regimes: 
Nonequilibrium equation-of-motion approach. 
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Inelastic effects in the Coulomb blockade and Kondo regimes of electron transport through molec- 
ular junctions are considered within a simple nonequilibrium equation-of-motion (EOM) approach. 
The scheme is self-consistent, and can qualitatively reproduce the main experimental observations 
of vibrational features in Coulomb blockade [H. Park et al., Nature 407, 57 (2000)] and Kondo 
[L. H. Yu et al., Phys. Rev. Lett. 93, 266802 (2004)] regimes. Considerations similar to the equi- 
librium EOM approach by Meir et al. [Phys. Rev. Lett. 66, 3048 (1991); ibid. 70, 2601 (1993)] are 
used on the Keldysh contour to account for the nonequilibrium nature of the junction, and dressing 
by appropriate Franck-Condon (FC) factors is used to account for vibrational features. Results of 
the equilibrium EOM scheme by Meir et al. are reproduced in the appropriate limit. 

PACS numbers: 73.23.Hk 72.10.Di 73.63.-b 85.65.-|-h 
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I. INTRODUCTION 

Fast development of experimental techniques in the 
area of molecular electronics makes it possible to ob- 
serve the response of molecular conduction junctions in a 
wide range of external parameters, such as source-drain 
and gate voltagesii Coulomb blockade (that characterizes 
the weak molecule- lead coupling limit), where transport 
through the molecular junction is suppressed due to high 
charging energy, and Kondo effect (encountered at suf- 
ficiently low temperature and strong molecule-lead cou- 
pling), when a correlation between localized (molecular) 
and band (contacts) electrons manifests itself in molec- 
ular junctions as a maximum in electrical conductance 
near Vsd ~ 0, were observed in the I/Vsd characteristics 
of such junctionsi^'^'^i^i^ii These are often accompanied 
by vibrational features that result from coupling between 
electronic and vibrational degrees of freedom. The latter 
can be associated with molecular center-of-mass motion^ 
or with intra molecular vibrations i^i^i^ 

Early theoretical approaches to transport in the 
Coulomb blockade regime were based either on linear 
response theory for near equilibrium situations^ i^^'^^i^^ 
or by treating transport at the level of quasi classi- 
cal rate equations F^^i^^ While the second approach to 
nonequilibrium transport is justified in the case of pure 
Coulomb blockade (where hopping between molecule 
and contacts is rare), the intermediate regime, e.g. 
the case of stronger molecule-leads coupling relevant 
for observation of nonequilibrium Kondo resonance, 
should be treated at a more sophisticated level. Re- 
cent approaches dealing with nonequilibrium Coulomb 
blockade and/or Kondo effect are based either on the 
slave-boson techniqu o^^i^^'^^i^^ , the equation-of-motion 
methodiiii^i^Si'^, the Fock-space rate equation schemoii, 
or the contour perturbation theory i-^-^i^^'^'^i^^i^^i^'^i^^'^^ In- 
elastic effects were not considered in the references above. 

Here we present a simple generalization of the 



equilibrium equation-of-motion approach used in the 
Coulombi2ii2£ regime (applied later also to the Kondo^ 
situation) to the case of nonequilibrium transport. 
The main difference between our approach and earlier 
nonequilibrium EOM studic o^^i^'^i^^ is a simple appealing 
structure of the Green function, the evaluation of which 
(in the absence of electron-phonon coupling) does not 
require a time-consuming self-consistent procedure. As 
was indicated earlier;^ this Green function expression re- 
duces to the exact solution both for an isolated molecule 
and in the limit of noninteracting electrons. We also 
generalize this basic scheme to include inelastic effects 
approximately, within an approach based on the Born- 
Oppenheimer approximation that is commonly used in 
Marcus theory of electron transferi^ Numerical calcu- 
lations are performed and qualitative correspondence to 
experimental data is demonstrated. 

Our model and theoretical procedure are presented in 
Section ini Numerical results for the Coulomb blockade 
regime are given and discussed in Section Hill The Kondo 
regime is discussed in Section ITVl Section IVl concludes. 



II. MODEL AND METHOD 



We describe the molecular junction within a single res- 
onant level (molecular electronic orbital) model, with 
electron-electron on-site repulsion (Hubbard term) and 
polaronic coupling to a local vibrational mode. The lat- 
ter is coupled to a bosonic thermal bath. The electronic 
orbital is coupled to two {L and R) free-electron reser- 
voirs representing the leads, each at its own equilibrium. 
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The corresponding Hamiltonian is 

K=L,RkeK,cr a- 
P K=L.RkeK,cr 

+ Uhini+ MQa E "'^ + E ^pQaQp 

where a =ti i is the electron spin index, Cka (Cfc^.) ^'"^ 
destruction (creation) operators for electronic state ka 
in the contacts, da (rfj-) destroys (creates) electron in the 
molecular orbital, d (a^) are second quantization opera- 
tors for the local vibrational mode, and bp (SJj) are the 
corresponding boson operators for thermal bath modes. 
Also 



1^ Qp^bf} + bl 



(2) 



are displacement operators for corresponding modes and 
fia ~ dlyda ■ Here and below we use h ~ \ and 
e = 1. After small polaron (canonical or Lang-Firsov) 
transformation^ the Hamiltonian takes the form (for de- 
tails see Ref. l3< 



K=L,RkeK,<j <y 

■T.^pbyb„+ J2 E (i4.cLrf".+H.c.) (3) 

/3 K=L.Rk£K,cr 
Unifii + ^ UpQaQp 



where 



and where 



U = U - 2MVwo 

Vkcr = VkcrXa 
Xa ^ CXp (iXaPa) ; K 



M 



(4) 
(5) 
(6) 

(7) 



is the phonon shift generator operator with 

Pa = -i{d- d^) (8) 

Pa, Eq.([8]), is the phonon momentum operator; we use 
the term phonon to characterize both molecular and bath 
vibrations. 

The Hamiltonian ([3]) is our starting point for the cal- 
culation of the steady-state current across the junction, 
using the nonequilibrium Green function (NEGF) expres- 
sion derived in Refs. |30||39| 



Here S^'^ are lesser/greater projections of the self- 
energy due to coupling to the contact K {K ^ L,R) 



^kAE) = ^fK{E)TKAE) (10) 
^kA^) ^ fKiEWKAE) (11) 

with ffciE) the Fermi distribution in the contact K and 

TkAE) = 27r ^ \Vka\^5iE - Sk) (12) 

keK 

The lesser and greater Green functions in ^ are 
Fourier transforms to energy space of projections onto 
the real time axis of the electron Green function on the 
Keldysh contour 

G.(ti,T2) = -^ < Tja{Tl)dl{T2) >H (13) 
= -^ < TjaiTi)Xa{n)dUT2)Xl{T2) >h 

where the subscripts H and H indicate which Hamilto- 
nian, ([T|) or ([31) respectively, determines evolution of the 
system, and Tc is the contour ordering operator. In what 
follows we use the second form and will drop the sub- 
script H while keeping in mind that time evolution is 
determined by the Hamiltonian We next decouple 
electron and phonon dynamics in the spirit of the Born- 
Oppenheimer theory within the Condon approximation 

Ga{Tl,T2) ^ GI^\ti,T2) K{ti,T2) (14) 

Gi'HT,,T2) = -l<Tja{Tl)dl{T2)> (15) 
K{tuT2) = <T,Xa{Ti)Xl{T2)> (16) 



where 



The shift generator correlation function K can be ex- 
pressed within the second order cumulant expansion in 
terms of the phonon Green function (for derivation see 
Ref. [37' 



K{ti,T2) = exp|A^ iDp^p^{Ti,T2)- < P^ > 

Dp^pAruT2) = -i< T,Pa{Ti)Pa{T2) > 



} (17) 
(18) 



while the phonon Green function D obeys approximately 
an equation which resembles the usual Dyson equation 



(19) 



^ dri ^ dT2 D^^Ip^ (t, n ) Hp„ P„ (n , T2 ) Dp^ P„ (r2 , r') 



with 



Hp„P„(ri,r2) = \Up?Dp,,p^,{T^,T2) 



(20) 
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(9) 



< E \yk.?[9kAr2,n)Gi'\n,T2)K{T^,T2) 

ke{L,R},(7 

+ (ti ^ T2)] 
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the analog of a self-energy, gk^a is the free electron GF 
in the contact, defined in (pij) below. 

To obtain an expression for the Green function G^f^ we 
follow the equation-of-motion (EOM) method of Meir, 
Wingreen, and Le o^°i'^° where it was applied for a near- 
equilibrium situation, except that we consider the EOMs 
on the Keldysh contour in order to take into account 
the nonequilibrium condition. In the spirit of the Born- 
Oppenheimer approximation we regard the shift genera- 
tor operators Xa as parameters incorporated into trans- 
fer matrix elements Vka- The solution of the electronic 
problem is thus carried out as in the absence of clectron- 
phonon couplin g^°'^° with renormalized parameters U 
{ U) and V { ^ V). The result is then averaged 
over the phonon subspace. This average is obviously 
not needed in the absence of electron-phonon coupling, 
M = 0, in which case G„ ^ G\^\ This leads to (for 
derivation see Appendix VK\ 



Gi^Hrur^) = [1- < >] G^^^ln, r2) (21) 



<n^>Gti{T,,T2) 



where the GFs cf] (« = {1, 2, 3, 4}) obey 



/ drG-i(ri,T)G(fj(r,r2) = <5(ti,T2) (22) 

J c 



with 



Gr.^(r,r') = 



(*|:-£^-t^) (23) 

-Y.^o{t,t') - J^asiT, t)] 

-S<,o(r,T') (24) 
+ U l^dT"Gtl{T,r")^^,{T",T') 



U l^ dT"G^:l{r,r")^,2iT",r') 



G-A{t,t') 



-I],o(T,T')-S.3(r,r')] 



(25) 



(26) 



given by 

S.o(r,r') =^ 114.1' <r,Xt(r)l„(T') > 
fc 

(27) 

SCTi(r, r') = < nfcg > 
k 

X [\Vk,f 9llir,r') <T,X,{r)XliT') > 

(28) 

+ \V,,\'gl^l{r,r') <TJlir)X,iT') > 
E,2(t, t') = S],3(t, t') - I],i(t, t') (29) 
S.3(r,r') = 9['l{r,r') < TMr)Xl{r') > 

k 

+ \V,,fg]^l{r,r')<T,Xlir)X,ir') > 

(30) 

with a denoting the spin opposite to cr. The free electron 
propagators g^^a and g^"*^, j = 1, 2 are defined by 



1-5 £ka 

OT 



OT 



^ ^ka ^cr 

OT 



5fe,.(T,T') = <5(r,T') 



(31) 



U 



&.r') = S{T,T') 



9i'iir,r')^6iT,T') (32) 



(33) 



For AI ~ O.V Franck-Condon (FC) factors (i.e. shift gen- 
erator correlation functions < X X^ > and < X'^ X >) 
should be taken as 1 in (gTll-dSni). Below the SEs in this 
case will be denoted Y^l^j (j = 0, 1,2,3). Note that the 
retarded projections of these are equivalent to the SEs 
introduced in Ref. [lol . For example, taking the retarded 
projection of (|30p and Fourier transforming to energy 
space leads to 



k£L,R 



1 



E + Eks — e<T — £s 
1 

E - £us - Ea + e 



u 



(34) 



Expressions for 'self-energies' S^-i {i = {0,1,2,3}) are 



which is identical to Eq.(9) in Ref. ^13 for i = 3. Other 
expressions are obtained in a similar way. 

Consider first the case with no electron-phonon cou- 
pling. The structure of expression (|2ip for the nonequi- 
librium GF g'^it^ is appealingly simple and has two impor- 
tant implications. First, it provides a convenient way for 
handling the Hubbard repulsion term U . While the case 
of weak electron-electron interaction can be handled by 
taking this term as a perturbation^ the case of strong 
interaction cannot be handled in this way, but includ- 
ing U in Hq makes standard diagrammatic techniques 
unusable 1^ This difficulty is circumvented by Ea. (pT|) . 
that expresses the system GF as a superposition (with 



4 



the level population n defining weight parameters) of 
simpler GFs associated with Hamiltonians that do not 
depend on U (apart from a parametric energy shift) for 
which the Wick's theorem is applicable. Secondly, by us- 
ing the EOM method on the Keldysh contour we are able 
to derive not only the retarded GF as in Ref. 1(1 but also 
the other projections, in particular the lesser GF that 
can be used to evaluate the level populations 



(35) 



This, together with Eg. ([21]) . lead to an explicit expression 
for < n„ >. Denoting 



one gets from ([2T|) 

< fla >= (1- < fla >)h,cr+ < n-ff > h,a 

and hence 

l-[h.S-h,s][h,a~hA 



< rirr > = 



(36) 



(37) 



(38) 



only qualitative description of the Kondo effect, since cor- 
relation between localized spin at the level and opposite 
spin cloud in the contacts is treated perturbatively. 

When electron-phonon interaction is present Eg. ([55]) 
remains valid. This results from the fact that K^{t,t) = 

1 so that G'^{t,t) ~ G^a^^ {t,t)] still, one has to deal 
with a self-consistent procedure. Indeed, the phonon 
GF Dp^p^ (and hence shift generator correlation func- 
tion /\ , see Eq. ((T7)) ) depends on the electronic GF G^f' 
through its 'self-energy' Ilp^p^, Eq. ((20l) . On the other 

hand, the electron GF Gir depends on the shift genera- 
tor correlation function K through its 'self-energies' J^^ri 
{i = {0, 1,2,3}), Eqs. (P7|) - ([5g|^ The resulting procedure 
is described in detail in Ref. [33. The only difference that 
enters here is the need to obtain the different self-energies 
defined in Eqs. (P7|) - pI)) . 

As discussed in Ref. [s^ the calculations involving 
electron-phonon interaction, when multiplication by the 
EG factor is necessary, are facilitated by repeatedly mov- 
ing between the time and energy domains. This is done 
using fast Fourier transform (EFT). In the calculations 
we use (followingi^) for the retarded projection of 'S,K,aO 



^(e) r 



.(0) ^(0) 



2s- £;x' 



(0) 



iW 



(0) 
K,cr 



(44) 



G-^j^ can be calculated from the Keldysh equation 

g[T{e) = Gif (£;) j:^<{e) G^t];{E) (39) 

with s|'2^ [i = 1,2,3,4) being lesser projections of the 

corresponding self-energies presented in Eqs.(l23l)-(l26l), We take w'^^^ = IOC/ and E^°\ taken at the Fermi level 



while its lesser projection is given by pop . where 



TkME) = -2Im S^:;o(^) 



(45) 



I.e. 



e1^) (r, r') = E.o(r, r') + S,3(t, t') (40) 
- S.o(r,r') -U f dT"G['liT,T")^^,ir",r') 

J C 



(41) 



4iir,r') - S.o(r,r') + U f dr" G[%,T")j:Mr" ,r') 

J c 

(42) 



defined to be the zero of energy {Ep = 0). This form will 
ensure convergence of the integrals. A band width ten 
times the Coulomb repulsion is enough to get essentially 
constant density of contacts states in the relevant energy 
region (wide band). F^''^ is taken much smaller than U 
to simulate the Goulomb blockade regime; exact numbers 
are indicated in calculation parameters below. 

The biased junction was characterized by the choice 



s5(r,r')-S.o(r,r') + S,3(r,T') 



(43) 



Ml 



V eVsd 



(1 - 77) eVsd (46) 



and expressions for So-i (« = {0,1,2,3}) given by (l??)) - 
dSO]). 

Since G^^j^ {i — 1,2,3,4) and therefore li^^ do not 
depend on < rio- >, Eg. ([55]) is an explicit expression for 
< n^r > and not, as might have expected, an equation 
that needs to be solved self-consistently. Eg. ((2T|) there- 
fore constitutes an explicit expression for Go- that can 
be evaluated directly once the G^-'^^ are known. Thus 
the Keldysh contour based consideration provides full 
information on the nonequilibrium system, and no sepa- 
rate considerations (as non-crossing approximation used 
in Ref. [sih are needed in order to estimate the level popu- 
lation. Note that both Ref. [Sll and our consideration give 



with voltage division factor 77 = 0.5. In calculations with 
M ^ 0, where an iterative procedure was used, con- 
vergence was assumed when population differences (elec- 
tronic population for both spins and vibrational popula- 
tion) between consecutive iteration steps were less than 
predefined tolerance, taken to be 10"**. The appHcation 
of a gate potential was represented by taking 



?^{Vg) ^ S^Vg = 0) + eVg 



(47) 



Note that Vg in (|T7)) is the effective potential at the 
molecule, which is usually considerably smaller than the 
bare potential applied to the gate. 

In what follows we apply the procedure outlined above 
in two situations. In section IIIII we focus on Coulomb 
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blockade phenomena. In section [TVl we describe the ap- 
phcation to Kondo physics by keeping the temperature 
low enough and by assigning finite lifetimes to the metal 
electrons. 



III. NUMERICAL RESULTS IN THE 
COULOMB BLOCKADE REGIME 

When dealing with the Coulomb blockade type calcu- 
lations, the electronic part (without the Franck-Condon 
(FC) factors) of the lesser and greater projections of So-j 
(j = 1,2,3) are obtained from Eas. (P51) - ([5(I)) and given 
by 

^i1''iE)= (48) 

K=L,R 

^T{E) = 

-i '^K,a{E^„)^K{E^a)\l- jK{Exa)\ (49) 



K=L,R 



+TKAE2a)lK{E2a)[l ~ fK{E2a) 



.(e)>,<, 



{E) = j:i'^>^<{E)-j:^:r^^{E) 



(e)>.<. 



(50) 



^.3 'iE) = 

I J2 [^KMEla)fK{Ei^)+rKAE2a)fK{E2a)] (51) 

^iTiE) = -^ E [^KAEi.)[l-fK{Ei.)] 

K=L,R 

+TkAE2.)[1- lK{E2a)]] (52) 

where Ei„ = Sa + Sg + U — E and E2a = E — -\- Sg. 
Retarded projection of the full SEs (after dressing by FC 
factors) are obtained using Lehmann representation!^ 

Consider first the situation where no electron-phonon 
coupling is present, M = 0. Figure [T}; shows a conduc- 
tance contour plot as a function of the gate and source- 
drain voltages for a system characterized by = —0.5, 
r^-*^ = 0.01, and T ~ 10^"* (all parameters are in units of 
U). Fig. [1^ presents average level population (soHd hue) 
and current (dotted line) plotted as a function of Vsd 
at fixed Vg = — [//4. I/Vgd curve shows two Coulomb 
addition plateaus, as is expected for a doubly degener- 
ate single level. Fig. [T|d is a similar graph as function 
of Vg at fixed Vsd = U /2. The usual Coulomb blockade 
diamond structure is observed in the bottom graph. Nat- 
urally, at high positive Vg the level is unpopulated, while 
at high negative Vg it is fully populated (< n >= 2). 
Within the conduction diamond the average population 
is 1, indicating the Coulomb blockade situation. Interme- 
diate regions provide fractional average populations due 
to partial occupation of the levels. 

The case Sa- e^, that may correspond to magnetic 
field removal of spin degeneracy is shown in Figure [51 We 
take the split levels to be = —0.6 and — —0.4, other 




015 




FIG. 1; (Color online) Elastic resonant tunneling. Average 
population (solid line, red; left axis) and current (dotted line, 
black; right axis) as function of (a) Vsd at fixed Vg = —U/4 
and (b) Vg at fixed Vsd = U/2. (c) Contour plot of dl/dVsd 
vs. Vg and Vsd- See text for parameters. Note that Vsd axis 
range in (a) goes beyond that in (c). 
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FIG. 2: (Color online) Elastic resonant tunneling under ap- 
plied magnetic field. Average populations of spin up (solid 
line, red) and spin down (dashed line, blue) levels vs. Vg at 
fixed Vsd = U/2 (top). Contour plot oidl/dVsd (bottom) vs. 
Vg and Vad- Inset shows two possible states (solid and dashed 
lines) of the molecule. See text for parameters. 

parameters are identical to those of Fig. [T] This split re- 
sults in splitting of the conductance lines as is shown in 
the bottom graph. Note the different intensity of the lines 
outside the diamond. The difference becomes even more 
drastic inside the diamond. This result is in agreement 
with experimental observatioui^ The calculated average 
population of the two spin levels (top graph), where again 
the source-drain voltage is fixed at Vsd = U/2, shows 
their complex dependence on gate voltage. This behav- 
ior can be understood within a simple argument. The 
molecule in the junction can be in either of the two states 
sketched in the inset of the top graph by solid and dashed 
lines. The observed average is the sum of the two con- 



tributions with weights representing probability for the 
system to be in the state. A qualitative explanation is 
based on the assumption that the system strives to be in 
a minimum energy situation (note that this explanation 
is only qualitative, since an energy minimum is not re- 
quired in the noncquilibrium transport case, however it 
might work to some extent in the blockade regime) . Thus 
the probability to be in the state indicated by solid lines 
in the inset is much higher than in the other. So, the most 
pronounced lines in conductance appear when chemical 
potentials cross the energy levels of this (solid line levels 
in the inset) state. Average population behavior can be 
explained with this consideration as well. 

In the presence of vibrational degrees of freedom inelas- 
tic co-tunneling (vibrational inelasticity) can be observed 
in conductanccj^ii!^ The situation is illustrated within a 
zero-order calculation^ using the parameters (in units 
oiU)T = 10-3, = -0.5, r^_^^ = 0.01, Loa = 0.2, and 
M = 0.4. The following points should be noted: 

1. Figurc[3^ shows the main Coulomb steps in the con- 
ductance map. In addition to elastic, vibrational 
sidebands corresponding to phonon creation by the 
tunneling electron are observed. Peaks correspond- 
ing to phonon absorption are not seen due to the 
low temperature employed in the calculation. 

2. Figure [3)3 represents the second derivative of cur- 
rent vs. source-drain voltage map. In addi- 
tion to resonant vibrational sidebands (lines along 
main Coulomb steps) observed in Fig. [3^ here one 
sees also inelastic electron tunneling spectroscopy 
(lETS) vibrational features (gate voltage indepen- 
dent off-resonant vibrational features) as well as 
weak lines corresponding to phonon annihilation. 

3. The absence of vibrational sidebands fro variable Vg 
for Vsd < tjJo is clearly seen from Fig. [3^. This issue 
was first addressed in Ref. [i^ and later confirmed 
by usi^ 

4. Suppression of the conduction signal at low 
source-drain voltage (the so called Franck-Condon 
blockade^) is seen from Fig. [3^ as well. At even 
stronger electron-phonon coupling (Figure (St; a 
zero-order calculation with the same parameters as 
in Fig. [3^ except that M = 0.6), the low voltage 
signal is suppressed completely. 

5. Note, that while experimentally the scales in Vg 
and Vsd where Coulomb blockade diamonds are ob- 
served are very different {Vsd is of order of Coulomb 
repulsion energy, 100 mV, while Vg spans 1 V), in 
our calculations they are comparable. The reason 
for this is that experimentally only part of the ap- 
plied gate voltage affects the position of the molec- 
ular level relative to contact Fermi energy. This 
is due to two reasons: first, capacitance factors 
(charging of the junction) play a role, and second, 
gate voltage can not be tuned to strongly affect the 
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molecule because of small sizes of the junctionj^ In 
our calculations however a rigid shift of molecular 
level is assumed. 



IV. THE KONDO REGIME 

The Kondo effecli^, a crossover from weak to strong 
coupling between localized (molecular) and band (con- 
tacts) electrons, manifests itself in molecular junctions as 
a maximum in electrical conductance near Vsd ~ at low 
temperatures. Conduction in this regime was described 
by Meir et ali^ within an EOM scheme. The treatment 
has focused on the retarded GFs, making it necessary 
to get level populations from a separate calculation us- 
ing the non-crossing approximation (NCA). In contrast, 
the NEGF EOM approach yields both the retarded and 
lesser GFs, and the needed level populations are obtained 
from the latter. This provides a single consistent theo- 
retical framework that, as we show below, reproduces 
the results of Ref. [3l|. It should be noted however that 
this approach is still an approximation, since truncating 
the EOM hierarchy, Eas. (|A12|) - (|A17p . implies neglect of 
correlations that may become important in the mixed va- 
lence situation when the level (shifted by Vg) is close 
to the Fermi energy. Therefore our nonequilibrium treat- 
ment of the Kondo regime is questionable beyond the low 
bias regime Ep — ^ eVsd, similar to the mean- field 
slave boson approac h^^i^^i^^i^^ where charge correlations 
are neglected by the mean-field approximation. In both 
approaches though the needed correlations in spin fluc- 
tuations are maintained; in the present approach this is 
done by keeping the correlation functions (jA6P - (jA8p as 
essential ingredients of the calculation. 

Consider first the purely electronic case, M = 0. 
Following^l we limit our consideration to the U —> oo 
limit. This leads to significant simplification while at the 
same time limiting the site to at most single occupancy 
as required for observation of the Kondo effect From 
Eqs. (ElD-dlSl) it follows that G^''^ - 1/C/ ^ in this 



limit, while G 
equation 



2,<T 



a 



(e,oo) 
2,cr 



satisfies the following Dyson 



i-^ ~ £a ] 5{ti,t) - Effo(ri,r) - t!'^\ti,t) 



(53) 



with SctO defined in ([27)) and from Eq. ([28|) (because 
g[}l ^ in the C/ ^ oo limit; c.f. Eq.l^) 

^l7\r,r')= J2 El^'=*l'(-'^-)3i'i(-'-') (54) 

K=L.B.k£K 

Thus from (HT)) it follows that the total GF in the [/ -> cx) 
limit is 

G(^'-) = [l-<n,>]G(:r^(Ti,r2) (55) 



FIG. 3: (Color online) Inelastic tunneling, (a) Contour plot 
of dl/dVsd vs. Vg and Vsd- (b) Contour plot of d?I /dV^d vs. 
Vg and Vsd- White regions correspond to values outside the 
scale; (c) Franck-Condon blockade. See text for parameters. 
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The Kondo peak diverges unless the finite Hfetime of 
metal electrons is taken into account. We incorporate 
this lifetime in the form introduced in Eq.(5) of Ref. [3l| 
(which associates lifetime with scattering off the molec- 
ular state). Note that the Lorentzian form adopted 
followingi^ for the coupling between molecule and con- 
tacts, Eq. (HH); prevents ultraviolet divergence of inte- 
grals such as (|B1[) and allows analytic evaluation of S^"^"* 
projections (see Appendix iBl Eqs. (jBlOp and (jBlip ). 
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FIG. 5: (Color online) d^I/dV^d for molecular junction in the 
Kondo regime. Shown are results for three choices of strength 
in coupling to vibration: M = r/2 (solid line, red), M = 2r/3 
(dashed line, blue), and M — T (dotted line, magenta). See 
text for other parameters. 



E(I) 



FIG. 4: (Color online) Bridge density of states in the Kondo 
regime for equilibrium (dashed line, blue) and nonequihbrium 
(solid line, red) situations. See text for parameters. 

Eqs. ([55|) - ([55|) lead to the following form for the re- 
tarded projection of 



G 



(e.oc) 

2,fT 



1- <nrr > 



(56) 



whereto (i^) and T.'-°^'{E) are defined in in Eqs. (gll) 
and (jBlOp respectively. These expressions arc identical 
to Eqs. (3) and(4) of Ref. '3^. Note however that < hg > 
is now calculated from the lesser projection 



Gi"'°°'<(S) 



[1- < fig. >] G. 



(e,oo)< 



[E) (57) 



Figure [3] presents the bridge density of states in equilib- 
rium (dashed line) and nonequihbrium (solid line) sit- 
uations. Parameters of the calculation arc (in units 



of r, 



(0) 



= r 



(0) 



^ R.a) 



T = 0.005, 



w^Z = 100 



S]- = ei = -z, 

As before the equilibrium Fermi energy 
defines the energy origin, and the noncquilibrium situa- 
tion is characterized by = Ep + \eV\ and ii^ ~ Ep. 
In equilibrium a Kondo peak at the Fermi energy is seen. 
It splits into two (at each of the electrode Fermi energies) 
when finite bias is applied. Comparing to Figs, la and b 
of Ref. [31I we see that the present formalism essentially 
reproduces these results. 

Inelastic effects are introduced into the picture as be- 
fore, by dressing transfer matrix elements by the shift 



operators, see Eqs. (HI]), and ((281). Figure [5] shows 
the result, obtained from such calculation for the sec- 
ond derivative of the current with respect to the source- 
drain voltage (bottom graph), for three choices of the 
electron-vibration coupling strength. Parameters of the 

calculation are (in units of F^P'') T — 0.025, £a = —2, 
= 100, Wo = 0.5. The solid, dashed, and dot- 
ted lines correspond to M = 0.5, 0.75, and 1, respec- 
tively. As is expected, increase in electron-vibration in- 
teraction destroys the Kondo effect. The reasons for this 
are (a) dephasing due to electron-vibration interaction 
and (b) shift of the energy level due to phonon reorga- 
nization. Electronic level shift downwards decreases the 
Kondo temperature {Tk ~ exp — 7r|ecr|/rCT see Ref. 15 ) 
thus destroying the Kondo peak. 

It should be emphasized that the vibrational struc- 
ture seen in Fig. O is a normal inelastic tunneling feature 
that is seen to persist also in the Kondo regime. This 
feature appears both in the Kondo and in the normal 
blockade regimes (see Figs. and [5]), as indeed was re- 
cently observed in the molecular junction experiment of 
Yu et ali^ The transition between these regimes (when a 
molecular orbital crosses the Fermi energy) can not be de- 
scribed by our approach for reasons outlined above. Also, 
Paaskc and Flensberg^ have recently applied a pcrturba- 
tive renormalization group to a limiting form of the same 
model in which the molecular electronic level is always in 
equilibrium with one side of the junction (the substrate 
in an STM configuration) and have shown that maintain- 
ing quantum coherence of vibrons, the effect disregarded 
in our treatment due to approximation p4p . may lead 
to enhancement of the exchange coupling and hence the 
Kondo temperature. 
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V. CONCLUSION 

We study inelastic effects in electron transport through 
a model molecular junction in Coulomb blockade and 
Kondo regimes. The approach is based on nonequi- 
librium generalization of the equation-of-motion scheme 
introduced by Meir et al.— >^ and is appealingly sim- 
ple. Inelastic effects are treated within a diabatic Born- 
Oppenheimer scheme. Important features of this ap- 
proach are correct analytical results for both isolated 
molecule (no contacts) and noninteracting {U = 0) cases, 
ability to reproduce results by Meir et al4i without ne- 
cessity of additional considerations to get the level pop- 
ulation, no necessity for self-consistency to get exact 
(within the scheme) results when the electron-vibration 
interaction is switched off, and unified treatment of both 
Coulomb and (to some extent) Kondo at noncquilibrium. 
The approach is able to reproduce experimental features 
qualitatively. 

Inelastic effects obtained within the model are reso- 
nant vibrational sidebands in the allowed, and lETS sig- 
nal in the blockaded, parts of the conductance map in 
Vg — Vsd coordinates, Franck-Condon blockade of trans- 
port for relatively strong electron- vibration interaction in 
the Coulomb blockade regime, and vibrational sidebands 
of the Kondo peak, as well as its quenching for strong 
vibronic coupling. 

Generalization of these considerations to the case of a 
two-site molecular bridge in the junction is straightfor- 
ward. The only problem is the large number of equations 
needed to be taken into account in this case. We post- 
pone such generalization for future study. 
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APPENDIX A: DERIVATION OF EQ. ((211) 

Here we derive Eq. (HI]). Note that the derivation 
does not depend on whether V ot V (and similarly U 
or U) is used for the system-leads coupling as long as the 
shift generator operator X is regarded as a scalar. We 
follow the procedure invented by Meir, Wingreen, and 
j^g Qio,30 fQj. ^Yie equilibrium situation and generalize it 
to the Keldysh contour case, in order to take into ac- 
count the noncquilibrium nature of molecular junction 
transport. During the derivation we will treat transfer 
matrix elements Vfccr, Eq. as numbers with the shift 
generator operators Xa, Eq. ([7]), incorporated into them 



as scalar parameters (a Born-Oppenheimer type approx- 
imation). However we'll have to keep track of their de- 
pendence on time (or more precisely contour variable) in 
order to get the phonon correlation functions K correctly 
at the end. 

We start from EOM for GF G^"^ (r, r'), Eq. USD, on 
the Keldysh contour 



.d_ 



G^^^ (r, r') = ^ (r, r') + UG^i^^^ (r, r') (Al) 



ke{L,R} 

new GFs on the r.h.s. have the form 

Ti'^ir,r') = -^<T,cUr)dUr')> 
G^^'\t,t') = -i< Tj4T)n4T)dUT') > 
Now we write EOMs for these GFs 



(A2) 
(A3) 



■ d 

^1=' 

OT 



^i'Jir,r') = VUr)Gi^\r,r') 



(A4) 

G^^'\T,T')^Sir,T') <n^> (A5) 
E [^L (-)ra!. (r, r') + {r)T^lll (r, r') 



£a-U 

OT 



While the EOM (|A4[) closes the chain of equations (its 

r.h.s. contains only G^f'), the EOM for G^^^-* yields new 
correlations in its r.h.s. defined by 



rSUr,r') = I < TAAr)n.ir)dUr') > 



(A6) 



^2l]Ar,r') = I < Tj,,{r)l{r)d,{r)dUr') > (A7) 
= ' < TccUr)dUr)dAr)dUr') > (A8) 
As a last step in the chain of EOMs we follow refer- 



encei 



.10.30 



by writing equations for the GFs (|A6P - (|A8 



i-^ £fco 

OT 



r[%{r,r') = V,4r)G(?^\ry) 



(A9) 



OT 



£a 



u 



vL{r)Gi"^Hr,r') (AlO) 

-E f^.MrSC(->-') + v^.Mrf^C(r,r') 



d 



i-^ £ks - £. 

OT 



£a 



V,,,{t) [G(^)(r,r')-G(2'=)(r,r')" 



(All) 



E [y^'^i-)^fS'kAr.T') - V,t.(r)r(5,,^(r,r') 
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On the right-hand-side of these equations we now have 
new, higher order GFs, T^'^^?^ ^ defined by the middle 
terms of Eqs. (Pl2)) - (PT7l) . GFs T^"^"^ and r^^^^' take ac- 
count of spin correlations in the leads. Closure of the 
(in principle infinite) EOM chain is achieved assuming 
that higher-order spin correlations in the leads can be 
neglected. Thus, following Ref. [13, the terms r'^*^^ are 
expressed in terms of lower order GFs 





= 


-I < Tj,,^{T)ck.{T)MT) 4 (r') >« 

(A12) 


p(3e) / 


r') = 


~t < TcCka{T)ck'a{T)dl{T) diir') >« 

(A13) 






-I < TA'a{T)cl{T)MT)diiT') >« 

(A14) 


p(3e) / 




(A15) 

Sk,k' <nka>G'-f{T,T') 


p(3e) / 




~l<T,Cks{T)cl,-^{T)d,{T)dl{T')> 

(A16) 


p(3e) ( 




-i < T,Cks{T)dUT)ck'a{T) diir') >w 

(A17) 



Now using (Pl2)) - (PT7)) in (|M)) - (PlT|) one can solve for 
rffcl (* = {1, 2, 3}) in terms of gS'^ and G^.^"' 



T 
T 
T 



1, fe,cr 
(2e) 

2, fc,(T 

(2e) 



(2e) 



< > Gi^) - Gi^-) 



(A18) 
(A19) 
(A20) 



where we have used short notation style with 'o' im- 
plying convolution of two functions on the contour [A o 
B){t,t')= J^dr" A{t,t")B{t",t'). These solutions are 

substituted into (jASp which gives G^^'^^ in terms of G^'^'' . 



Finally, the last result together with (jA4p can be used in 
(|Aip to get equation for G^*^' in the form 



Gi' 



G 



(e) 



U<h,> Gii o GZ 



(A21) 



G\';l {i = {1,2,3,4}) are defined in Eqs. ([22ll-(|26ll. 
while 'self-energies' entering these definitions are given 
by Eqs. dllll-dSOl). 

In order to simplify the structure we rewrite it in the 
form 



G'-^) = [1- < n, >] Gii+ <n,> {g^^ + UG[% o G[i] 

(A22) 

and note that 
{...] = G^^ioG^i\G-,l 



U 



(A23) 



The last equation follows from Gj^ ^Gj ^ = G4^G3^. 
Substitution of (|M3)) into (|M2)) leads to jHI). The re- 
tarded projection of (pijl is the final result of Ref. [l3- 



APPENDIX B: ANALYTICAL EXPRESSION FOR 
SELF-ENERGY S^'^' 



Here we derive analytical expressions for retarded and 
lesser projections of S^^'', Eq. ((54)) . under Lorentzian 
assumption for coupling between molecule and contacts, 
Eq. (j44|l . In the case of a dense continuum of states in the 
contacts (assumed here) the sum in ([M)) can be converted 
to an integral, then retarded and lesser projection of the 
SE (in energy domain) are 



K=L.R' 



s(-)<(£;) 



E - e - e„ + £s + i^isl^ 
(Bl) 



+ 00 



K=L.,R-^-°° 2^ (£; - e - + e^)' + (7^/2)^ 
de rK,s(e)/if(e) 



E 

<=L,1 

^ E 

K=L,R 



where second line of (|B2p is correct for the case of T ^ 
(relevant for observation of the Kondo peak) . 
Introducing 



(B2) 



x = I3{t- (Ik) 
we arrive at integrals of the form 

1 1 



dx ■ 



-00 



dx ■ 



{x - xi){x - X2){x ~ X3) + 1 
1 



1 



(a; — xi){x — X2){x — X3){x — X4) + 1 
for (|B1[) and (jB2p respectively, where 

' AO) , ,w(o) 



(B3) 

(B4) 
(B5) 



Xi 






X2 






X3 




'e 


X4 


— X3 





(B6) 
(B7) 

y) (B8) 
(B9) 



with K ~ L, R. These integrals can be taken analyt- 
ically by complex contour integration, the poles are at 
xi, X2, X3 {x4 in the case of integral (jBSp ). and also at 
Un = «7r(27i + 1): n = 0, ±1, ±2, . . . Performing the in- 
tegration one arrives at the following expressions for the 
SE projections 
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K=L,R 



47r 



^2 + + 1,12) E2 - - 7*/2), 



/c.cr Ai.cr 



(BIO) 



p(0) 






27r 


[^2-*(<i-7./2) 







4 i;2 + *(<i + 7*/2) 



(0) 



K=L,R 



27r 



Im 



(Bll) 



(0) 



-Im 



^2-*(<i-7./2) ^;2+^(<i+7*/2) 



W, 



(0) 



(0) , 75 



2^^l + (<i+75/2P, 



with E2 = E — Ea+eg — E^i^l and where ip is a Psi (digamma) function.— Note that it is the second term in Eq. (jBlOp 
which is responsible for Kondo effect appearance. ^ 
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